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The joint NASA-ESA mission LISA relies crucially on the stability of the three spacecraft constel- 
lation. All three spacecraft are on heliocentric and weakly eccentric orbits forming a stable triangle. 
It has been shown that for certain spacecraft orbits, the arms keep constant distances to the first 
order in the eccentricities. However, exact orbitography shows the so-called 'breathing modes' of 
the arms where the arms slowly change their lengths over the time-scale of a year. In this paper 
we analyse the breathing modes (the flexing of the arms) with the help of the geodesic deviation 
equations to octupole order which are shown to be equivalent to higher order Clohessy- Wiltshire 
equations. We show that the flexing of the arms of LISA as given by the 'exact' solution of Keplerian 
orbits, which gives constant armlengths to the first order in eccentricity and whose maximum flex- 
ing amplitude is ~ 115, 000 km, can be improved, by tilting the plane of the LISA triangle slightly 
from the proposed orientation of 60° with the ecliptic to obtain a maximum flexing amplitude of 
~ 48, 000 km, reducing it by a factor of ~ 2.4. The reduction factor is even larger if we consider 
the corresponding Doppler shifts, for which the reduction factor reaches almost a factor of 6. We 
solve the second order equations and obtain the general solution. We then use the general solution 
to establish the optimality of the solutions that we have found. 



I. INTRODUCTION 



Detection of low frequency (below 100 mHz) gravitational waves (LFGW) is a challenging endeavour because 
terrestrial instruments are dominated by seismic noise in this frequency range. There is not much hope of isolating 
the noise in this frequency range, which mainly arises from the direct Newtonian coupling of ground motions with 
the test masses. This is why space missions have been given serious considerations from early on, because LFGW are 
extremely important in the physics of and around massive black holes which are hosted by most of galaxies, including 
ours. The LISA mission [l| involves three spacecraft forming a triangle of side ~ 5 x 10 6 km, linked by optical 
beams. The gravitational wave (GW) signal is read from the beat note generated between an incoming beam and 
the local laser, so that six such signals are generated between the three pairs of spacecraft. In a classical Michelson 
interferometer, the requirements on the laser phase noise are strongly reduced by the symmetry between the two arms. 
In the case of LISA, a comparison of the laser beam from a distant spacecraft with the local oscillator would demand 
a relative frequency stability better than the GW amplitude, which is quite unrealistic. This was realised from the 
very beginning. The idea which allows relaxing the requirement of frequency stability is to digitally reproduce an 
interferometric configuration by mixing the six elementary data flows with delays corresponding to virtual optical 
paths between spacecrafts. The corresponding technique is called Time-Delay Interfcrometry (TDI) But the 
practical realisation of a triangular configuration orbiting the sun leads to a structure only approximately rigid. In 
fact, the three spacecraft are on almost circular orbits of small eccentricity e, and the mutual distances are constant 
during the year only to the first order in e. When an exact calculation of the mutual distances between pairs of 
spacecraft is performed using Keplerian orbital equations, a 'breathing' motion of the arms appears, normally termed 
as 'flexing' of the arms, with amplitude of about 115,000 km. This relative motion between spacecraft causes a Doppler 
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effect, but at frequency (~ 1CP 7 Hz) out of the detection band 10~ 4 — 10 _1 Hz. It has nevertheless been shown Q 
that this motion makes the delay operators of TDI, variable in time and thus prevents the exact cancellation of the 
laser frequency noise. One can however, suppress the noise but at the cost of complex application of TDI techniques 
involving noncommuting delay operators. When first generation TDI combinations are used, the residual noise is 
proportional to the time derivative of the flexing, so that a significant reduction of the slow orbital Doppler effect 
would make it necessary to reconsider the TDI strategy. We have found it useful, therefore, to discuss the stability of 
the LISA triangle to higher orders in e instead of only to the first order, in order to investigate whether better orbits 
are possible, with reduced breathing amplitude. 

The paper is organised as follows. In section |n] we firstly show that a small correction to the angle between the 
LISA plane and the ecliptic is able to reduce by a factor ~ 2.4 the RMS flexing amplitude and more importantly the 
peak to peak Doppler shift by a factor ~ 5.6 which is relevant for TDI. In section ITTT1 we show how a generalisation 
of the Clohessy- Wiltshire equations up to second order, namely, the octupolar expansion of the Newtonian potential 
and the related geodesic deviation, allows us to reproduce almost exactly the behaviour as obtained from the exact 
orbits. This paves for us a natural way for introducing extra degrees of freedom in the orbital model in order to 
establish the optimality of the solution we have found. In section llVl we achieve this by obtaining the general solution 
to the second order equations and vary the arbitrary constants to arrive at the optimal results. In this analysis, for 
simplicity, we have ignored the effects of the gravitational pulls of Jupiter and Earth which are of the same order. We 
note that the flexing occurs mainly due to two reasons: (i) the non-existence of exact Keplerian orbits of spacecraft 
allowing constant arm-lengths, and (ii) the perturbation due to planets in our Solar system. In this paper we address 
the first issue only. 



II. REDUCING THE FLEXING OF THE 'EXACT' SOLUTION BY SLIGHTLY TILTING THE PLANE 

The LISA mission requires that the distances between the three spacecraft remain constant to the first order in the 
parameter a — 1/2R. a is proportional to the eccentricity e to the first order: a = V3e to this order in e. We take 
I = 5 x 10 6 km and R = 1 A. U. ~ 1.5 x 10 8 km, and thus a ~ 1 /60. The constancy of armlengths alone does not lead to 
a unique choice of the orbits. It has been shown that all the points around a circular reference orbit-in a plane making 
an angle of ±-| with respect to the plane of reference orbit satisfies the above condition^, 0, @, Q, 0, However, 
all these orbits have different arm length variation to the second order in a, also known as flexing. In reference 0] 
(henceforth referred to as paper I), we have presented a set of elliptical orbits satisfying the above requirement. Here 
we extend our analysis to second order in a which describes the flexing of the arms quite accurately. We show this 
by comparing the second order results with the exact results which agree to within few parts in 10 4 . Thus the second 
order analysis allows us to handle the flexing in a convenient and reliable manner. The eventual aim of this excercise 
is to minimise the flexing which could be important for the exploitation of LISA. 

For this purpose we need to decide on the measure of flexing. There are several ways of characterising this quantity. 
We can take this to be the peak to peak variation of the armlength - the difference between the maximum distance 
between two given spacecraft and the minimum distance between the two spacecraft over the period of one year; or 
we can define it as a r.m.s variation of armlength over the period of a year. An important measure is the Doppler 
shift which depends on the relative velocity between two spacecraft. We investigate all these measures of flexing in 
this paper. We also make the simplifying assumption that the orbits of the spacecrafts 2 and 3 are obtained by rigidly 
rotating the orbit of spacecraft 1 by 120° and 240° (this is described more precisely in the next subsection). Due to the 
symmetry assumed in this model, the measures of flexing although different from each other are independent of the 
spacecraft pair we choose. Therefore without loss of generality we consider spacecraft 1 and spacecraft 2. If Z12 (i) is 
the distance between the spacecraft 1 and 2 at time t, then the peak to peak flexing is equal to maxZi2(£) — imnl^it), 



FIG. 1: The figure shows the geometry of the orbits and of LISA. The barycentric frame is labelled by {X, Y, Z} while the CW 
frame is labelled by {x,y,z}. SCI, SC2 and SC3 denote the three spacecraft. The radius of the reference orbit is taken to be 
R = 1 A. U. and S denotes the Sun. 



where the maxima and minima are taken over the period of one year. The ^23 (^) and li3(t) are time-shifted versions 
of li2(t) by one third and two thirds of an year respectively. The goal of this section is to explore the possibility of 
reducing the flexing by tweaking the orbital parameters of the spacecraft orbits. One way is to vary the tilt of the 
plane of LISA and study the effect on the arm length variation as a function of the tilt angle and then minimise 
the flexing by appropriately choosing the tilt. We show that the flexing can be reduced by appropriately choosing 
the tilt and the flexing can be reduced to a minimum. The generic problem involves exploring the space of orbital 
parameters of the three spacecraft orbits. Because of the symmetry assumed here, we are restricting ourselves to a 
subspace, whose dimensionality is smaller by a factor of three than the dimensionality of the full parameter space of 
orbits. Following paper I, we start with the exact orbital equations. 



The exact orbits 



We choose the barycentric frame with coordinates {X, Y, Z} as follows: The ecliptic plane is the X — Y plane and 
we consider a circular reference orbit of radius 1 A. U. centred at the Sun. Let 5 be a small correction of the order 
a, to the tilt of the LISA plane. We choose the axes so that the tilt is 7r/3 + 5. We choose spacecraft 1 to be at its 
highest point (maximum Z) at t — 0. This means that at this point, t — and Y = 0. The orbit of the first spacecraft 
is an ellipse with inclination angle e, eccentricity e and satisfying the above initial condition. The geometry of the 
configuration is shown in Figure ^ From the geometry, e and e are obtained as functions of 5, 
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and the orbit equations for the spacecraft 1 are given by: 
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(1) 



X\ = i?(cosi/'i + e) cose, 
Yi = R\J 1 — e 2 sin -01, 
Zi = R(costpi + e) sine. 



(2) 
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FIG. 2: (a) Plot of distance £12 as a function of time and tilt angle S. (b) Plot of the optimal variation in arm length as 
function of time. For comparison we have also given I12 for tilt equal to ^. 



The eccentric anomaly ipi is implicitly given in terms of t by, 

ipi + e sin tpi = Sit , (3) 

where t is the time and O is the average angular velocity. The orbits of the spacecraft 2 and 3 are obtained by rotating 
the orbit of spacecraft 1 by 27r/3 and 47r/3 about the Z— axis; the phases V'2 J '/'3 ) however, must be adjusted so that 
the spacecraft are about the distance I from each other. The orbital equations of spacecraft k = 2, 3 are: 

Xk = X\ cos <7fc — Y\ sin a k , 
Yfe = X\ sin <7fc + Y\ cos a k , 

Z k = (4) 
where Ofc = (fc — 1) with the caveat that the "01 is replaced by the phases tpk, where they are implicitly given by, 

tpk + esin^ fe = fit - (k - T)— = Clt - 07.. (5) 

These are the exact orbital equations for the three spacecraft. 

We investigate the arm length variation as function of the tilt angle 5. The distance between spacecraft 1 and 2, 
I12 = {(Xi - X 2 ) 2 + (Yi - Y 2 ) 2 + (Z x - Z 2 ) 2 } 1/2 , as a function of t and 6 is shown in Figure 2 (a). Figure 2 (b) shows 
constant 8 sections of the 3 D plot in (a). For (5 = the peak to peak variation in armlength is about 115,000 km 
while the r.m.s. variation is about 36000 km. However, as can also be seen from the plots, the variation in armlength 
can be reduced by a factor little more than 2. Our goal is to determine the tilt angle which minimises the flexing 
and also to determine this minimum. Our second goal is to show that sufficiently accurate results can be obtained by 
going only to the second order in a. We start from the exact orbits and then make a second order expansion in a in 
the next subsection. 
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The orbits to second order in a 



In this section, we obtain the expression for orbits to the second order in a. From the exact expressions for the 
inclination e and the eccentricity e given in Eq. flj. to second order in a we obtain, 



a a 2 



_ aS , 

a 2 aS . . 

e = a -7i + 7r (6) 

Eq. JJJJ can be solved iteratively to the required order in e to obtain an expression for tpk(t) as a function of time, 
substituting for e from Eq. © and retaining terms only to the second order in a. We then have, 

2 

o> a 

i>k{t) = 4>k - sin^fc + — (sin20 fe - 3sin<?!> fc ) + a<5 sin0 fc , (7) 

where 4>k = fit — ^r(fc — 1) = fit — <7fc, for the fcth spacecraft. The orbital equation for the fcth spacecraft can be 
obtained by substituting these values in Eq. Q . The expressions take on simpler form if we transform to frames tied 
to LISA. Here we use these approximate, albeit, simpler equations in order to compute the inter-spacecraft distances. 
The first transformation is to the Clohessy- Wiltshire (CW) frame. 

Clohessy and Wiltshire 01 make a transformation to a frame - the CW frame {x, y, z} which has its origin on the 
reference orbit and also rotates with angular velocity fi. The x direction is normal and coplanar with the reference 
orbit, the y direction is tangential and comoving, and the z direction is chosen orthogonal to the orbital plane. They 
write down the linearised dynamical equations for test-particles in the neighbourhood of a reference particle (such as 
the Earth). Since the frame is noninertial, Coriolis and centrifugal forces appear in addition to the tidal forces. 

We take the reference particle to be orbiting in a circle of radius R with constant angular velocity fi. Then the 
transformation to the CW frame {x, y, z} from the barycentric frame {X, Y, Z} is given by, 

x = (X - R cos fit) cos fit + (Y - R sin Qt) sin Vlt , 
y = - (X - R cos fit) sin fit + (Y - R sin fit) cos fit , 

z = Z. (8) 

The second transformation to the LISA frame is given by rotating the CW frame about y axis by an angle ir/3 + 6. 
With these two transformations the orbital equations for the fcth spacecraft to the second order in a are given by, 

2 1 
x' k = —=Ru cos^fc + — i?a 2 [13 — 6cos(/)fe — 7cos2<?!)fe] , 
v3 24 

2 1 

y' k = -=Ra sin^fc + 2Ra6 sin 0& + —Ra [sin20fc — 3sin^fc] , 

V 3 3 

4 = -^Ra 2 [11 - 10 cos0 fe -cos2^ fc ] . (9) 
8%/ 3 

Using these equations it is easier to compute the distance between spacecraft. 



C. Optimal tilt angle 

Here we obtain the expression for the distance between the spacecraft as function of tilt angle 5. This can be 
computed in a straight forward way from Eq. (1) - (5) in the barycentric frame. The computations require less effort 
in the CW frame or in the {x 1 , y', z'} frame. In the {&', y', z'} frame it is easy to see that the z' contribution in the 
inter-spacecraft distance drops out since it is of higher order in a. As remarked earlier we just need to consider Z12, 



6 




0.4 0.5 O.f 

Time in years 



FIG. 3: The variation of the lengths of the arms of LISA (the breathing modes) is shown in units of 10 6 km. The dashed curve 
shows the variation in I12 over a year computed to the second order in a, while the solid curve shows the exact variation. The 
discrepancy between the curves is less than 0.03 %. 



the distance between SC 1 and SC 2, because of the symmetry assumed in the model. To the second order in a we 
get the following expression for Z12 : 



h 2 (t,5)=l + Al 12 (t,6) 



(10) 



where, 



AZ12M) 



a 2 R 
16\/3 



48 1 ? 



S 1 



15 cos 6> + 48 



5 A 



cos 29 — cos 36* 



(11) 



Here 5 = aSi and 9 = fit — tt/'S. Since we restrict 5 ~ o(a), Si ~ o(l). This expression is quite simple and easy 
to manipulate in order to extract the relevant information. Again, we see from figure 0] the exact orbits (solid curve) 
and the second order (dashed curve) approximation in a match very well. This fact motivates the analytic derivation 
of the equations of motion to second order in a. We do this in section HTT1 We make the following observations: 

• To the first order in a, l\i is constant in time and equal to I. This result is known in earlier literature and was 
also obtained in paper I from the CW equations. It is only from the second order in a that the variation in 
arm-length appears. 

• The variation in arm-length obtained up to the second order in a is close to the exact variation in arm-length. 
First we compare the variations in arm- lengths as a function of time for the tilt angle of ^ or when 5 — 0. When 
S = we get from Eq. Ijll[). 



Al 12 (t,5 = 0) = ^-al U -5 cos 9 + 10 cos 261- icos36» 



(12) 



With a — 1/60 we have plotted li2(t, 6 = 0), depicted by the dashed curve, in Fig. |3| The dominant term is of 
double frequency - the term in 2flt - which gives two cycles in one year. The heights of the peaks are modulated 
by the single frequency term in fit. Also the curve is 'pushed' up by the constant term or the DC term. The 
term in 3flt is small and has little effect on the overall qualitative behaviour. The solid curve shows Z12 (t, (5 = 0) 
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for the exact Keplerian orbits. It is important to note that the second order result is very close to one obtained 
from the exact Keplerian orbits. The maximum difference in the flexing is about 0.03 percent. 

• We plot the peak-to-peak amplitude of the flexing over the year as a function of S. This is shown in the Figure 
El As is seen here, there is a range of values of S for which the flexing is optimal, namely, 0.0083 < S < 0.0125. 
In this interval the curve appears flat; we show that the curve is flat in this region when computed to the second 
order in a and the exact model is numerically found to be similar. This can be inferred from Figure 2 (b) where 
we have plotted l\2 as a function of t for various values of S, particularly when S lies in the optimal range. From 
Figure 01 it is observed that the difference between maximum and minimum I12 remains constant ~ 48000 km, 
in the optimal range of i5. This result can be obtained exactly from Eq. (jll[) which is correct to second order in 
a as follows: 

Writing x = cos 9 and expanding the cos 29 and cos 39 in terms of cos 9, we obtain, 

a 2 R ' 



AZ 12 



4^3 



x 3 - 24 - - S 1 )x 2 + 3x + 3 



(13) 



The extrema of AZ12 with respect to time can be obtained by setting the derivative of A/12 with respect to x 
equal to zero. This yields an equation for x parametrized by 5±: 

x 2 - 2(5 - 8Si)x + 1 = . (14) 

By evaluating the discriminant of this equation, we find that it has real roots only if S\ lies out of the interval 
(0.5,0.75). Within this interval the roots are complex and the extrema are obtained by setting x = ±1; it is 
clear from Eq. (|13|) that the minimum of I12 (or AZ12) is attained when x = 1 and the maximum when x = — 1. 
The peak-to-peak amplitude is just the difference between these two extrema and is easily computed from Eq. 
(|13fl . Wc find from these considerations that the peak to peak amplitude in the optimal region is: 

2a 2 R al , % 

AZi 2 , max - AZ 12 , min = — j=- = -j= ~ 48, 000km . (15) 
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We also compute the r.m.s. variation in the arm-length from Eq. i|ll|) as follows: 



(A/ 12 ) = V3a 2 R ( - - S 1 ) , 



(16) 



where the bracket denotes the r.m.s. variation over 1 year. Clearly only the constant term contributes to the 
average. The variance of l 12 denoted by V(li2) is easily obtained by taking the sum of the squares of the 
coefficients of the time dependent terms because of the mutual orthogonality of the terms, and then dividing by 
2. The division by the factor of 2 is because (cos 2 k9) — 1/2, k — 1, 2, 3. Thus, 



n/i 2 ((5 1 )) = (/ 12 -(/ 12 )) 2 (,5 1 ) = 



a 4 R 2 
1536 



225 + 2304 - - <5i 



1 



(17) 



Clearly the variance is minimum when Si =5/8 which gives 5 = aSi = 0.0104 which is right in the centre of the 
optimal interval [0.5,0.75] of S\. From this value of Si = 5/8, the minimum variance as well as the minimum 
r.m.s. amplitude of the breathing mode can be computed. The minimum r.m.s. amplitude of flexing is ~ 16, 000 
km. The gain factor g in the r.m.s. amplitude with respect to the case S — is easily obtained from Eq. l(T7|l 
to be, 



fJ 

With Si = 0, the peak to peak amplitude is: 

A^12,max — A/ 12min = 
whereas in the optimal range for Si it is reduced to: 

.max ,min 



' ^12(31=0)) 

V(h 2 (Si = 5/8)) 



2.23. 



(18) 



, i \/3 r 

!x-x — 
R 2 



4V6 - 9 



115,000 km, 



; 1 1 

/ x — x 



48, 000 km . 



(19) 



(20) 



R 2^3 

which is a factor of about 2.4. The quantities we try to minimize are however the Doppler shifts, i.e. the time 
derivatives of the arm-length variations. The time derivative of A^i 2 is: 



d Al - 
~dt 12 16^3 



15sin6»-96 - - Si) sin 29 + 3 sin 39 



We can also compute the variance of d/dt(Ali2): 



V 



dt 



Al u 



a 4 R 2 fl 2 
1536 



234 + 9216 [--Si 

8 



(21) 



(22) 



It is obviously a minimum for Si = 5/8, but we note that the gain in Doppler shift is now: 



9d 



3834 
234 



4.05, 



(23) 



so that it is even larger than the gain on the flexing itself. If we now investigate the peak to peak value of the Doppler 
shift, we get about 45.8 m/s for the standard (S = 0) case, and only 8.2 m/s for the optimised case, so that the ratio 
is about 5.6, which seems interesting to note. The relative velocity is shown in Figure where its peak to peak value 
is apparent. It is seen from the figure that there are long periods of constant Doppler shifts for the optimised case. 
These results could be useful in the context of TDI. 
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FIG. 5: The relative velocity ^(^12) is shown for the two cases 5i = (dashed curve) and the optimised case 5i = 5/8 (solid 
curve) over the course of a year. For the optimised case the gain in peak to peak variation in relative velocity is ~ 5.6. 

III. THE EQUATIONS OF MOTION OF SPACECRAFT TO OCTUPOLE ORDER 

A. Geodesic deviation equation to second order in the separation vector 

We begin with the Newtonian geodesic deviation equation to the first order in the separation vector: 

fj o + Jf 6 V = 0, (24) 

where r/ a (t), are the components of the 3-vector separating the two particles at any time t and K% — The 
$(r) = — k/r, k = GM, is the Newtonian gravitational potential, where r is radial coordinate with the centre of force 
at the origin, G the Newtonian gravitational constant and M the mass of the Sun. Here the commas denote partial 
derivatives with respect to the coordinates. Upper and lower indices have been used to facilitate the summation 
convention. 

If the reference orbit is circular with radius R, the geodesic deviation equation in the barycentric frame is given by: 



ij° + lL(R*6S-3X a X b )ri b = 0, 



(25) 



where X a is the vector (X, Y, Z). In vector notation (a bar over the quantity denotes a vector quantity and a hat 
denotes a unit vector), 



r/ + fl 2 f}-3fl 2 R{R.fj) =0, 



(26) 



where fl 2 = k/R 3 and R denotes a unit vector in the radial direction. 

As seen before, the CW frame is a suitable frame. In this frame, the separation vector fj is related to the CW 
coordinates by, 



x = rf cos Q,t + rf sin fit , 
y = —rf sin fit + r\ v cos fit , 

z = rf, 



(27) 
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V° + + ^'lv b r) c = 0, (29) 



®'l = SK^bXc + S a c X b + S bc X a ) + —X a X b X c . (30) 



Then transforming Ea. H26(l to CW frame, it is found to be identical to the CW equations: 

x - 2ny - m 2 x = , 
y - 2m = , 

z + n 2 z = 0. (28) 

However, in order to obtain the flexing of the arms we need to go to higher order in a. We go to the second order 
in a or which is equivalent to the second order in separation vector fj. The second order equation contains octupolar 
terms of the gravitational field which essentially drive the flexing of the arms. The Taylor expansion of the Newtonian 
gravitational potential to the third order leads to the following equation for the separation vector fj: 

1 

— < 

2 

where 

^X c + S a c X b + S bc X a ) + ^ 

Or in the vector form, the above equation can be written as: 

- - 3Q 2 3 15 

fj + Q 2 fj - 3n 2 R(R.fj) - — (R.fj)fj - + ^n 2 R(R.fj) 2 = 0. (31) 

The octupole terms in the equation are quadratic in the separation vector. These equations take on a simple form 
when transformed to the CW frame described earlier in Eq. (JBJ. 
In the CW frame, Eq. H31[) simplifies to: 

302 

£-2n y ~3n 2 x + — {2x 2 -y 2 -z 2 ) = 0, 
2R 

m 2 

y + 2flx —xy = 0, 

S + fl 2 z xz = 0. (32) 

R 

These equations can also be thought of as higher (second) order CW equations. 

In order to display the order of the various terms in a, it is useful to employ dimensionless coordinates. We define 
the dimensionless variables in units of / as, 

x ~ y ~ z ~ p 

x=-,y = -,z = -,p=-, (33) 

where p 2 = x 2 + y 2 + z 2 . 

We also set f2 to be unity for convenience. In the dimensionless variables, the equations (|32|l reduce to: 

x - 2y - 3x + 3a{3x 2 - p 2 ) = , 
y + 2x — 6axy = , 
z + z — 6axz = . (34) 

B. The Equations of Motion via the Lagrangian 

The equations of motion can also derived with the Lagrangian for a central force expanded around a point on the 
reference orbit. The Lagrangian is: 



L 2 



\ (X 2 + Y 2 + Y 2 ) + * . (35) 

2 V ) (X 2 + Y 2 + Z 2 ) 1/2 
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Here the mass of the test particle is ignored; it is not important for the central force problem. However if control 
terms and force terms are added then the mass of the test particle must be included. Here we ignore these effects. 
We carry out our analysis in the CW frame. The Lagrangian given by Ea. (|35[) in the CW- frame takes the form: 



L = 



i {i 2 + f + z 2 + (R + x) 2 + y 2 } + y (R + x) 



xy 



(R + x) +y 2 + z< 



1/2 



(36) 



Then the force free Lagrangian equations of motion are given by: 



x-2y- (1- (3){R + x) = 0, 
y + 2x - (1 - 0)y = 0, 
z + (3z = 0, 



(37) 



where, 



2x ( X \ 2 



(I) 



(- 

\R 



(38) 



Expanded to the octupole order, the quantity (3 is given by: 



a < 3a; 1 
^ 1 -R + 2 



3 <1 



31 * 



(39) 



To this order of approximation, Eq. I|37|) are identical to Eq. I|32|) . 



IV. THE SOLUTIONS TO THE EQUATIONS OF MOTION 

We first obtain a general solution to the equations of motion in terms of arbitrary constants and then describe the 
family of solutions with tilt 5 ~ o(a), which matches with the exact Keplerian orbits Eq. J5J) to second order in a. 
Finally, we show that the variation in arm-length is optimal in two ways: (i) r.m.s sense and (ii) peak to peak sense. 
This is achieved by optimally choosing the arbitrary constants in the general solution. These results are valid in the 
neighbourhood of the exact orbits considered in section III Al 



A. The General Perturbative Solution around 5 — 

Accordingly, we write, x = xq + axi, y = yo + ayi, z = zq + az±, where, Jo, yo and Sq are solutions to the equations 
of motion to quadrupole order. These solutions satisfy certain criteria which we describe below. 

To the first order in a the xo,yo,zo satisfy the CW equations J2HJ- As shown in paper I, the solutions of these 
equations satisfying no drift, no offset and which have the test-particle equidistant from the origin are the following: 

x = ip cos(t + 0o ) , 
y = -posing + <f> ) , 

z = — Po cos(< + (/>o ) , (40) 

where po and 0o are fixed for a specific orbit. For the LISA spacecraft, p~o = -^=. The orbits of the spacecraft are 
easily obtained by the appropriate choice of O . For convenience in solving the equations, we choose 0o to be zero. 
The orbits with non-zero 0o are easily obtained by translation. 
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To the second order in a, we get the equations for Xi,yi and z\ as, 

xi - 2yi - 3xi + 3(3x1 — pg) = , 
yi + 2xi - 6x yo = , 

Z\ + Z\ — 6io^o = 0. (41) 

We first solve the y-equation. We substitute the solutions for xq and yo from the CW solutions so that the 

y-equation becomes: 

3 

yi + 2xi + -pi sin 2t = . (42) 

Integrating we have, 

3 

yi + 2xi - -pi cos 2t = A, (43) 

where A is an integration constant. We substitute yi in the equation for xi to obtain: 

1 5 

xi + xi = - cos 2 t + 2A + - . (44) 



Then integrating we get, 



1 5 

£i = B cost + Csini - 777 cos 2t + 2A+ - . (45) 
24 8 



Finally, integrating the equation for Z\ ,we get, 



1 \f3 

Z\ = Ecost + Fsint ^cos2tH . (46) 

4%/3 4 

After integrating the yi-equation, namely, Eq. 143(1 , we collect the solutions: 

1 5 

xi = B cos t + Csint cos2t + 2A + - , 

24 8 

1 5 

yi = -2Bsint + 2Ccost + -sin2t- 3At t + D, 

6 4 

1 \/3 

zi = Ecost + Fsint ^cos2£H . (47) 

4V3 4 

Here, A, B, C, D, E and F are the integration constants and (|47() is the general solution for the perturbations to 
second order in a. These constants can be determined from initial conditions. However, the term linear in t in the 
solution for yi represents an unbounded drift. Since such solutions imply instability, the drift must be removed. This 
is achieved if we set A — — j|. 



B. The family of solutions with varying tilt near 5 = 

In this subsection we show that the family of solutions with varying tilt S ~ o(a) is consistent with the general 
solution obtained in Eq. (|47|l . For this purpose we consider the exact orbits of the spacecraft given by Eq. - Q- 
We are therefore choosing a specific set of orbits for which the spacecraft distances are constant and equal to I to the 
first order in the eccentricity (or equivalently a). If these orbits are expanded to second order in a (note 5 = aSi 
where 5i ~ o(l)), the following equations can be deduced for spacecraft 1 in the CW frame in which the general 
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solution Eq. (|47() was obtained: 



Ra 2 
x = -^cost + 2Ra 

y = -2-^ sint + 2 J Ra 2 
v3 



1 _ 5i 
4 2 



Si 



1 „ 5 
C0St_ 24 C ° S 24 



1 



sin t H — sin 2t 
6 



22a cos t + 2Ra 



1 /, 1 „ V3 

— = (oi — llcosf = cos2t H — — 

2^ 4^3 4 



(48) 



Thev match with Eg . (|4Tjl if the constants are chosen as follows: A = — yj, 5 = j — |<Ji, C = 0, D = 0, 25 = (^ 1_ 1) 
and F = 0. 

Constructing the full solution (adding first and second order solutions) from Eq. H40|) and Eq. (|47H . we obtain 
identical solutions, 



I 



I 2 



r^COSt , 

2V3 2i? 



1 Si 



4 2 



1 



— — — cos t — — cos It — — 



24 



24 



5i I sint H — sin 2t 

" ' 6 



I I 2 
Zl = 2 C ° St+ 2R 



1 ,r 1 „ \/3 

— p= (Oi — 1) cost ^ cos2t H — — 

2V3 4V3 4 



(49) 



where we have now switched back to coordinates having dimensions of length and x±, y±, z\ being the coordinates of 
spacecraft 1. The coordinates of spacecraft 2 and 3 in the CW frame are obtained by replacing t by t — 2tt/3 and 
t — Att/3 respectively. 



C. Establishing the optimality of the solution 



We have shown that the second order CW equations are sufficient for an accurate description of the breathing of 
LISA. We are therefore allowed to use the general solution found in section llV Al for the purpose of optimisation. If 
we compute the distance between SC 1 and SC 2 using the general solution, with the same arbitrary constants for all 
SC, we find: 



1X2 = I 



1+^(5(55 + ^2?) -A COS*- 



-{3B-V3E) cos 29 



-i(3C-V3F)sin2#- — 
2 V ' 48 



The time varying part of ?x2 is now: 



«X2 



(I 



12 



cos 36 



16 



i (32? - V3E) cos 29+l(3C- ViF) sin 29 
2 2 48 



(50) 



(51) 



Since the terms are orthogonal, the minimisation of the variance amounts to minimising each coefficient separately. 
The arbitrary constants, namely, B, C, E, F occur only in the coefficients of the cos 28 and sin 28 terms while the rest 
of the terms have constant coefficients. Thus the minimisation amounts to setting the coefficients of cos 28 and sin 28 
equal to zero. This gives the minimum r.m.s. variation of about 16,000 km which was the same result obtained at 
the end of section [D] 

If we now introduce different constants B,C,E,F for each SC, it is easy to show that only terms in cos k6, k = 1,2,3 
and sin 28 appear in the distance, in such a way that only the coefficients of cos 28 and sin 28 can be made to vanish. 
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The result is then exactly as the preceding one, showing that this is indeed a true r.m.s. optimum in the neighborhood 
of the exact solution that we have assumed. This is a general result. 

We also minimise the peak to peak flexing by considering the general solution. First, we consider the particular 
case when 3C — \/3F = 0. Then as in section ITU 

(52) 

where x = cos 6. Following the analysis analogous to section [TJ] we set the derivative of AZ12 with respect to x equal 
to zero which yields the quadratic equation: 

x 2 - 8(35 - V~3E)x + 1 = . (53) 

This equation has complex roots when —0.25 < 35 — ^/3E < 0.25 implying that the extrema of AZ12 are attained 
when x = ±1 and when this condition is satisfied. From Eq. I|52l) one obtains the minimum peak to peak variation 
as a l/VZ ~ 48000 km identical to the one found in section [ij] However, this family is larger, having two parameters 
5 and E than the one considered with varying tilt which has just one parameter 5. In particular, we verify that for 
the family of solutions with varying tilt 5 investigated in section HVBI the condition —0.25 < 35 — y/3E < 0.25, after 
substituting their values in terms of Si, leads to, 0.5 < 81 < 0.75 in agreement with the result obtained earlier in 
section ITTI 

If we do not apriori assume that 3C - V3F is zero, then the simple analysis as performed above, fails. We therefore 
resort to numerical methods. We first observe that the time varying part of I12 up to the constant Za\/3/4 from Eq. 
ijSTJ is given by: 

5 1 

F(a,b;9)=a cos 29 + b sin 29 - - cos 9 - — cos 36* , (54) 

8 24 

where, a = 35 — V3F and b = 3C — \/3F. The peak to peak value of F we denote by AF. We observe that F 
satisfies the following symmetries: F(—a, b; 9) = —F(a, b; it — 8) and F(a, —b; 8) = F{a, b; —9) implying AF(a, —b) = 
AF(—a, b) = AF(a, b). Thus AF is an even function of a, b. We can therefore restrict ourselves to non-negative values 
of a and b when searching for the optimum value of AF. Also for large a, b, AF ~ 2y 'a 2 + b 2 , so that the minimum 
of AF must lie in the neighborhood of the origin with a,b ~ o(l). A straight forward numerical computation shows 
that the minimum value of AF occurs when —0.25 < a < 0.25 and b = 0, giving the minimum peak variation of 
~ 48, 000 km as before. This is again a general result. 



AZ 12 = 



24 



12(35 - V3E)x 2 +3x- 125 



V. CONCLUDING REMARKS 

In this work we have studied the variation of the lengths of the arms of LISA in the neighbourhood of the exact 
solution that gives constant armlcngths to first order in a. We investigate the flexing of arms of LISA by going to 
the second order in the parameter a. We numerically establish that the second order results provide a near accurate 
description of the orbits and the flexing of the arms. Because the second order expressions are tractable and amenable 
to analytic treatment, a lot of information can be gleaned from the expressions themselves. We show that by slightly 
changing the tilt angle of the plane of LISA from the proposed one of 60° with the ecliptic, one can reduce the flexing 
of the arms. This is investigated with the aid of the exact Keplerian orbits as well as by analysing to second order 
in a. We then obtain the tilt angle for which the flexing of the arms is minimum. We find that the peak to peak 
amplitude of flexing is reduced by a factor of about 2.4. More importantly, in the context of application of TDI 
techniques, we show that the peak to peak Doppler shift is reduced by a factor of 5.6. 

Secondly, we obtain the equations of motion to second order in a from the Newtonian geodesic deviation equation 
to the quadratic order in the separation vector, which involves the octupolar expansion of the Newtonian gravitational 
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potential of the Sun. These equations can also be obtained with a Lagrangian approach. The general solutions to 
these equations are found perturbatively, in the neighbourhood of the first order solution for which the armlcngths are 
constant to first order in a. We then check that the family of solutions with varying tilt is reproduced from the general 
solution by choosing the arbitrary constants appropriately. Finally we use the general solution to establish that our 
solution of the minimum flexing of the arms is more generally valid in the neighborhood of the exact solution assumed. 
This is important from the point of view of a realistic launch because the flexibility in the choice of parameters allows 
for some error, which is inevitable, and this is possible without sacrificing performance. 



Acknowledgments 



SK would like to acknowledge DST, India for the WOS-A. 



[1] "LISA: A Cornerstone Mission for the Observation of Gravitational Waves", System and Technology Study Report, 2000. 
[2] M. Tinto and J. W. Armstrong, Phys. Rev. D 59, 102003 (1999); M. Tinto and S. V. Dhurandhar, Living Reviews in 

Relativity, http://www.livingreviews.org/lrr - 2005 - 4, (gr-qc/0409034). 
[3] M. Tinto, F. B. Eastabrook and J. W. Armstrong, Phys. Rev. D 69, 082001 (2004). 
[4] S. V. Dhurandhar, K. R. Nayak, S. Koshti and J-Y. Vinet, Class. Quantum Grav., 22, 481 (2005). 

[5] M.A. Vincent and P.L. Bender, Proc. Astrodynamics Specialist Conference (Kalispell USA) (Univelt, San Diego, 1987) 
1346. 

[6] W. M. Folkncr, F. Hcchlcr, T. H. Sweetser, M. A. Vincent and P. L. Bender, Class. Quant. Grav., 14, 1405 (1997). 

[7] S. P. Hughes, Preliminary optimal orbit design for Laser Interferometer Space Antenna, 25th Annual A AS Guidance and 

Control Conference, Breckenridge CO, Feb. 2002. 
[8] F. Hechler and W. M. Folkner, Adv. Space Res., 32, 1277 (2003). 
[9] T.H. Sweetser, Class. Quantum Grav., 22, S 429 (2005). 
[10] W. H. Clohessy and R. S. Wiltshire, Journal of Aerospace Sciences, 653 - 658 (1960); 

D. A. Vallado, Foundations of Astrodynamics and Applications, 2nd edition 2001, Microcosm Press Kluwer; 

also in S. Nerem, http://ccar. Colorado. edu/asen5050/lecturel2.pdf(2003). 



Reducing the flexing of the arms of LISA 



K. Rajesh Nayak 1 , S. Koshti 2 , S. V. Dhurandhar 2 and J-Y. Vinet 3 
The University of Texas at Brownsville, 
80 Fort Brown, Brownsville, TX 78520, USA. 
2 IUCAA, Postbag 4i Ganeshkind, Pune - J^ll 007, India. 
3 ARTEMIS, Observatoire de la Cote d'Azur, 
BP 4229, 06304 Nice, France. 
(Dated: July 24, 2005) 

The joint NAS A-ES A mission LISA relies crucially on the stability of the three spacecraft constel- 
lation. All three spacecraft are on heliocentric and weakly eccentric orbits forming a stable triangle. 
It has been shown that for certain spacecraft orbits, the arms keep constant distances to the first 
order in the eccentricities. However, exact orbitography shows the so-called 'breathing modes' of 
the arms where the arms slowly change their lengths over the time-scale of a year. In this paper 
we analyse the breathing modes (the flexing of the arms) with the help of the geodesic deviation 
equations to octupole order which are shown to be equivalent to higher order Clohessy- Wiltshire 
equations. We show that the flexing of the arms of LISA as given by the 'exact' solution of Keplerian 
orbits, which gives constant armlengths to the first order in eccentricity and whose maximum flex- 
ing amplitude is ~ 115, 000 km, can be improved, by tilting the plane of the LISA triangle slightly 
from the proposed orientation of 60° with the ecliptic to obtain a maximum flexing amplitude of 
~ 48, 000 km, reducing it by a factor of ~ 2.4. The reduction factor is even larger if we consider 
the corresponding Doppler shifts, for which the reduction factor reaches almost a factor of 6. We 
solve the second order equations and obtain the general solution. We then use the general solution 
to establish the optimality of the solutions that we have found. 



I. INTRODUCTION 



Detection of low frequency (below 100 mHz) gravitational waves (LFGW) is a challenging endeavour because 
terrestrial instruments are dominated by seismic noise in this frequency range. There is not much hope of isolating 
the noise in this frequency range, which mainly arises from the direct Newtonian coupling of ground motions with 
the test masses. This is why space missions have been given serious considerations from early on, because LFGW are 
extremely important in thephysics of and around massive black holes which are hosted by most of galaxies, including 
ours. The LISA mission [l| involves three spacecraft forming a triangle of side ~ 5 x 10 6 km, linked by optical 
beams. The gravitational wave (GW) signal is read from the beat note generated between an incoming beam and 
the local laser, so that six such signals are generated between the three pairs of spacecraft. In a classical Michelson 
interferometer, the requirements on the laser phase noise are strongly reduced by the symmetry between the two arms. 
In the case of LISA, a comparison of the laser beam from a distant spacecraft with the local oscillator would demand 
a relative frequency stability better than the GW amplitude, which is quite unrealistic. This was realised from the 
very beginning. The idea which allows relaxing the requirement of frequency stability is to digitally reproduce an 
interferometric configuration by mixing the six elementary data flows with delays corresponding to virtual optical 
paths between spacecrafts. The corresponding technique is called Time-Delay Interferometry (TDI) Q. But the 
practical realisation of a triangular configuration orbiting the sun leads to a structure only approximately rigid. In 
fact, the three spacecraft are on almost circular orbits of small eccentricity e, and the mutual distances are constant 
during the year only to the first order in e. When an exact calculation of the mutual distances between pairs of 
spacecraft is performed using Keplerian orbital equations, a 'breathing' motion of the arms appears, normally termed 
as 'flexing' of the arms, with amplitude of about 115,000 km. This relative motion between spacecraft causes a Doppler 



effect, but at frequency (~ 10 -7 Hz) out of the detection band 1CT 4 — 10 _1 Hz. It has nevertheless been shown |3| 
that this motion makes the delay operators of TDI, variable in time and thus prevents the exact cancellation of the 
laser frequency noise. One can however, suppress the noise but at the cost of complex application of TDI techniques 
involving noncommuting delay operators. When first generation TDI combinations are used, the residual noise is 
proportional to the time derivative of the flexing, so that a significant reduction of the slow orbital Doppler effect 
would make it necessary to reconsider the TDI strategy. We have found it useful, therefore, to discuss the stability of 
the LISA triangle to higher orders in e instead of only to the first order, in order to investigate whether better orbits 
are possible, with reduced breathing amplitude. 

The paper is organised as follows. In section |nj we firstly show that a small correction to the angle between the 
LISA plane and the ecliptic is able to reduce by a factor ~ 2.4 the RMS flexing amplitude and more importantly the 
peak to peak Doppler shift by a factor ~ 5.6 which is relevant for TDI. In section III we show how a generalisation 
of the Clohessy- Wiltshire equations up to second order, namely, the octupolar expansion of the Newtonian potential 
and the related geodesic deviation, allows us to reproduce almost exactly the behaviour as obtained from the exact 
orbits. This paves for us a natural way for introducing extra degrees of freedom in the orbital model in order to 
establish the optimality of the solution we have found. In section IV, we achieve this by obtaining the general solution 
to the second order equations and vary the arbitrary constants to arrive at the optimal results. In this analysis, for 
simplicity, we have ignored the effects of the gravitational pulls of Jupiter and Earth which are of the same order. We 
note that the flexing occurs mainly due to two reasons: (i) the non-existence of exact Keplerian orbits of spacecraft 
allowing constant arm-lengths, and (ii) the perturbation due to planets in our Solar system. In this paper we address 
the first issue only. 



II, REDUCING THE FLEXING OF THE 'EXACT' SOLUTION BY SLIGHTLY TILTING THE PLANE 

The LISA mission requires that the distances between the three spacecraft remain constant to the first order in the 
parameter a = 1/2R. a is proportional to the eccentricity e to the first order: a — y/3e to this order in e. We take 
I = 5x 10 6 km and R — 1 A. U. ~ 1.5 x 10 8 km, and thus a ~ 1/60. The constancy of armlengths alone does not lead to 
a unique choice of the orbits. It has been shown that all the points around a circular reference orbitjin a plane making 
an angle of ±f with respect to the plane of reference orbit satisfies the above conditional Q, 0, Q, fl Q. However, 
all these orbits have different arm length variation to the second order in a, also known as flexing. In reference |J| 
(henceforth referred to as paper I), we have presented a set of elliptical orbits satisfying the above requirement. Here 
we extend our analysis to second order in a which describes the flexing of the arms quite accurately. We show this 
by comparing the second order results with the exact results which agree to within few parts in 10 4 . Thus the second 
order analysis allows us to handle the flexing in a convenient and reliable manner. The eventual aim of this excercise 
is to minimise the flexing which could be important for the exploitation of LISA. 

For this purpose we need to decide on the measure of flexing. There are several ways of characterising this quantity. 
We can take this to be the peak to peak variation of the armlength - the difference between the maximum distance 
between two given spacecraft and the minimum distance between the two spacecraft over the period of one year; or 
we can define it as a r.m.s variation of armlength over the period of a year. An important measure is the Doppler 
shift which depends on the relative velocity between two spacecraft. We investigate all these measures of flexing in 
this paper. We also make the simplifying assumption that the orbits of the spacecrafts 2 and 3 are obtained by rigidly 
rotating the orbit of spacecraft 1 by 120° and 240° (this is described more precisely in the next subsection). Due to the 
symmetry assumed in this model, the measures of flexing although different from each other are independent of the 
spacecraft pair we choose. Therefore without loss of generality we consider spacecraft 1 and spacecraft 2. If Zi2(i) is 
the distance between the spacecraft 1 and 2 at time t, then the peak to peak flexing is equal to maxZi2(t) — minimi), 



FIG. 1: The figure shows the geometry of the orbits and of LISA. The barycentric frame is labelled by {X, Y, Z} while the CW 
frame is labelled by {x,y,z}. SCI, SC2 and SC3 denote the three spacecraft. The radius of the reference orbit is taken to be 
R = 1 A. U. and S denotes the Sun. 



where the maxima and minima are taken over the period of one year. The ^23 (*) and 1 13(f) are time-shifted versions 
of /12(f) by one third and two thirds of an year respectively. The goal of this section is to explore the possibility of 
reducing the flexing by tweaking the orbital parameters of the spacecraft orbits. One way is to vary the tilt of the 
plane of LISA and study the effect on the arm length variation as a function of the tilt angle and then minimise 
the flexing by appropriately choosing the tilt. We show that the flexing can be reduced by appropriately choosing 
the tilt and the flexing can be reduced to a minimum. The generic problem involves exploring the space of orbital 
parameters of the three spacecraft orbits. Because of the symmetry assumed here, we are restricting ourselves to a 
subspace, whose dimensionality is smaller by a factor of three than the dimensionality of the full parameter space of 
orbits. Following paper I, we start with the exact orbital equations. 



A. The exact orbits 

We choose the barycentric frame with coordinates {X, Y, Z} as follows: The ecliptic plane is the X — Y plane and 
we consider a circular reference orbit of radius 1 A. U. centred at the Sun. Let S be a small correction of the order 
a, to the tilt of the LISA plane. We choose the axes so that the tilt is n/3 + S. We choose spacecraft 1 to be at its 
highest point (maximum Z) at t = 0. This means that at this point, t = and Y — 0. The orbit of the first spacecraft 
is an ellipse with inclination angle e, eccentricity e and satisfying the above initial condition. The geometry of the 
configuration is shown in Figure From the geometry, e and e are obtained as functions of 5, 
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and the orbit equations for the spacecraft 1 are given by: 



1/2 



(1) 



X\ = R(costpi + e) cose, 



Y\ = Ry 1 — e 2 sin^i, 
Z\ = R(costpi + e) sine. 



(2) 



(a) 



Time in years 
(b) 



FIG. 2: (a) Plot of distance £12 as a function of time and tilt angle S. (b) Plot of the optimal variation in arm length as 
function of time. For comparison we have also given I12 for tilt equal to -| . 

The eccentric anomaly ipi is implicitly given in terms of t by, 

tpl + e sin ipi = fit , (3) 

where t is the time and O is the average angular velocity. The orbits of the spacecraft 2 and 3 are obtained by rotating 
the orbit of spacecraft 1 by 2n/3 and 47r/3 about the Z— axis; the phases *p2,ip3, however, must be adjusted so that 
the spacecraft are about the distance I from each other. The orbital equations of spacecraft k = 2, 3 are: 

X k = Xi cos (T fe - Yi sin a k , 
Yfe = Xi sin (T fe + Y\ cos a k , 

Z k = Z x , (4) 
where a k = (k — 1) with the caveat that the tpi is replaced by the phases V'fe, where they are implicitly given by, 

ip k + es'mipk = fit - (k - 1)— = fit - a k . (5) 

These are the exact orbital equations for the three spacecraft. 

We investigate the arm length variation as function of the tilt angle S. The distance between spacecraft 1 and 2, 
I12 = {{Xi - X 2 f + (Yi - Y 2 f + (Zi - Z 2 ) 2 } 1/2 , as a function of t and S is shown in Figure 2 (a). Figure 2 (b) shows 
constant S sections of the 3 D plot in (a). For S = the peak to peak variation in armlength is about 115,000 km 
while the r.m.s. variation is about 36000 km. However, as can also be seen from the plots, the variation in armlength 
can be reduced by a factor little more than 2. Our goal is to determine the tilt angle which minimises the flexing 
and also to determine this minimum. Our second goal is to show that sufficiently accurate results can be obtained by 
going only to the second order in a. We start from the exact orbits and then make a second order expansion in a in 
the next subsection. 



B. The orbits to second order in a 



In this section, we obtain the expression for orbits to the second order in a. From the exact expressions for the 
inclination e and the eccentricity e given in Eq. to second order in a we obtain, 

a a 2 . 
a 2 aS 

Eq. (2J can be solved iteratively to the required order in e to obtain an expression for tpk{t) as a function of time, 
substituting for e from Eq. © and retaining terms only to the second order in a. We then have, 

9 

.a a 
4>k(t) = 4>k - sin^fc + — (sin2</> fc - 3sin(f> k ) + a5 sin0 fc , (7) 

where cf>k = Qt — =£-(k — 1) = Qt — <7fc, for the fcth spacecraft. The orbital equation for the fcth spacecraft can be 
obtained by substituting these values in Eq.(@J. The expressions take on simpler form if we transform to frames tied 
to LISA. Here we use these approximate, albiet, simpler equations in order to compute the inter-spacecraft distances. 
The first transformation is to the Clohessy- Wiltshire (CW) frame. 

Clohessy and Wiltshire ^jj] make a transformation to a frame - the CW frame {x, y, z} which has its origin on the 
reference orbit and also rotates with angular velocity Q. The x direction is normal and coplanar with the reference 
orbit, the y direction is tangential and comoving, and the z direction is chosen orthogonal to the orbital plane. They 
write down the linearised dynamical equations for test-particles in the neighbourhood of a reference particle (such as 
the Earth). Since the frame is noninertial, Coriolis and centrifugal forces appear in addition to the tidal forces. 

We take the reference particle to be orbiting in a circle of radius R with constant angular velocity Q. Then the 
transformation to the CW frame {x, y, z} from the barycentric frame {X, Y, Z} is given by, 

x = (X - R cos Qt) cos fit + (Y ~ R sin fit) sin Qt , 
y = - (X - R cos Sit) sin Qt + (Y - R sin Qt) cos Qt , 

z = Z. (8) 

The second transformation to the LISA frame is given by rotating the CW frame about y axis by an angle 7r/3 + 5. 
With these two transformations the orbital equations for the fcth spacecraft to the second order in a are given by, 

2 1 
x' k = —=Ra cos(f>k + — i?a 2 [13 — 6cos(/)fe — 7cos2<?!)fe] , 
v3 24 

2 1 

y' k = -=Ra sin^fc +2RaSsm<f>k + izRa 2 [sin20fc — 3sin^fc] , 

v3 3 

4 = ~^Ra 2 [11 - 10 cos0 fc -cos2</> fe ] . (9) 
8%/ 3 

Using these equations it is easier to compute the distance between spacecraft. 



C. Optimal tilt angle 

Here we obtain the expression for the distance between the spacecraft as function of tilt angle 6. This can be 
computed in a straight forward way from Eq. (1) - (5) in the barycentric frame. The computations require less effort 
in the CW frame or in the {a/, y', z'} frame. In the {x 1 , y', z'} frame it is easy to see that the z' contribution in the 
inter-spacecraft distance drops out since it is of higher order in a. As remarked earlier we just need to consider Zi 2 , 



5.1 




Time in years 

FIG. 3: The variation of the lengths of the arms of LISA (the breathing modes) is shown in units of 10 6 km. The dashed curve 
shows the variation in £12 over a year computed to the second order in a, while the solid curve shows the exact variation. The 
discrepancy between the curves is less than 0.03 %. 



the distance between SC 1 and SC 2, because of the symmetry assumed in the model. To the second order in a we 
get the following expression for l i2 : 



h 2 {t,S) = l + Al 12 (t,5) 



(10) 



where, 



16\/3 



48 - -£i - 15cos 



5 A 



cos 29 — cos 30 



(11) 



Here 5 — a5\ and 6 = Qt — tt/3. Since we restrict 5 ~ o(a), 5± ~ o(l). This expression is quite simple and easy 
to manipulate in order to extract the relevant information. Again, we see from figure 0] the exact orbits (solid curve) 
and the second order (dashed curve) approximation in a match very well. This fact motivates the analytic derivation 
of the equations of motion to second order in a. We do this in section III. We make the following observations: 

• To the first order in a, Z i2 is constant in time and equal to I. This result is known in earlier literature and was 
also obtained in paper I from the CW equations. It is only from the second order in a that the variation in 
arm-length appears. 

• The variation in arm-length obtained upto the second order in a is close to the exact variation in arm-length. 
First we compare the variations in arm- lengths as a function of time for the tilt angle of ? or when 5 — 0. When 



5 = we get from Eq. (fTTf) - 



V3 



Al 12 (t,6 = 0) = -i—al I 6-5cos6> 



10cos2f3- -cos36> 



(12) 



With a — 1/60 we have plotted li 2 (t, 5 — 0), depicted by the dashed curve, in Fig. The dominant term is of 
double frequency - the term in 2£lt - which gives two cycles in one year. The heights of the peaks are modulated 
by the single frequency term in fit. Also the curve is 'pushed' up by the constant term or the DC term. The 
term in 3Sli is small and has little effect on the overall qualitative behaviour. The solid curve shows h 2 (t, 5 — 0) 




for the exact Keplerian orbits. It is important to note that the second order result is very close to one obtained 
from the exact Keplerian orbits. The maximum difference in the flexing is about 0.03 percent. 

• We plot the peak-to-peak amplitude of the flexing over the year as a function of S. This is shown in the Figure 
Q] As is seen here, there is a range of values of S for which the flexing is optimal, namely, 0.0083 < S < 0.0125. 
In this interval the curve appears flat; we show that the curve is flat in this region when computed to the second 
order in a and the exact model is numerically found to be similar. This can be inferred from Figure 2 (b) where 
we have plotted Z12 as a function of t for various values of 6, particularly when S lies in the optimal range. From 
Figure it is observed that the difference between maximum and minimum Z12 remains constant ~ 48000 km, 
in the optimal range of 6, This result can be obtained exactly from Eq. (1 1 H which is correct to second order in 
a as follows: 

Writing x = cos 9 and expanding the cos 29 and cos W in terms of cos 9, we obtain, 



a 2 R 
4V3 



x 3 - 24 ( - - (Si ) x 2 + 3x + 3 



(13) 



The extrema of A/i 2 with respect to time can be obtained by setting the derivative of AZi 2 with respect to x 
equal to zero. This yields an equation for x parametrized by Si: 



„2 



X 



2(5-8<5i)x + l = 0. (14) 



By evaluating the discriminant of this equation, we find that it has real roots only if 8\ lies out of the interval 
(0.5,0.75). Within this interval the roots are complex and the extrema are obtained by setting x = ±1; it is 
clear from Eq. 113t that the minimum of l\% (or AZ12) is attained when x = 1 and the maximum when x = — 1. 
The peak-to-peak amplitude is just the difference between these two extrema and is easily computed from Eq. 
Ill ■'■St . We find from these considerations that the peak to peak amplitude in the optimal region is: 

2a 2 R al 

A/i 2 , max - AZ 12 , min = = — ~ 48, 000km . (15) 



We also compute the r.m.s. variation in the arm-length from Eq. itTTjl as follows: 



(A/ 12 ) = V3a 2 R H-SA , 



(16) 



where the bracket denotes the r.m.s. variation over 1 year. Clearly only the constant term contributes to the 
average. The variance of h 2 denoted by V(li2) is easily obtained by taking the sum of the squares of the 
coefficients of the time dependent terms because of the mutual orthogonality of the terms, and then dividing by 
2. The division by the factor of 2 is because (cos 2 k9) = 1/2, k = 1, 2, 3. Thus, 



V(l 12 (6 1 )) = (l 12 -(l 12 )) 2 (5 1 ) = 



a 4 R 2 
1536 



225 + 2304 



Si 



1 



(17) 



Clearly the variance is minimum when Si =5/8 which gives 8 = aSi = 0.0104 which is right in the centre of the 
optimal interval [0.5,0.75] of 8%, From this value of <$i = 5/8, the minimum variance as well as the minimum 
r.m.s. amplitude of the breathing mode can be computed. The minimum r.m.s. amplitude of flexing is ~ 16, 000 
km. The gain factor g in the r.m.s. amplitude with respect to the case 8 — is easily obtained from Eq. I|17J1 
to be, 



9 = 

With Si — 0, the peak to peak amplitude is: 

Aii2,max — AZi2 im m = 

whereas in the optimal range for Si it is reduced to: 

.max ,min 



' V(h2(Si = 0)) 
V(h 2 (8i = 5/8)) 



2.23. 



(18) 



l X — X 

R 2 



4V6- 9 



115,000 km, 



/ 1 

x — x 



R 2y/3 



48, 000 km . 



(19) 



(20) 



which is a factor of about 2.4. The quantities we try to minimize are however the Doppler shifts, i.e. the time 
derivatives of the arm-length variations. The time derivative of Ali 2 is: 



— Al - 

«ft 12 Tevf 



15 sin 9 - 96 



Si sin 20 + 3 sin 30 



We can also compute the variance of d/dt(Ali 2 ): 

V 



1536 



234 + 9216 ( - - 8t 



(21) 



(22) 



It is obviously a minimum for i5i = 5/8, but we note that the gain in Doppler shift is now: 



9d 



3834 
234 



4.05, 



(23) 



so that it is even larger than the gain on the flexing itself. If we now investigate the peak to peak value of the Doppler 
shift, we get about 45.8 m/s for the standard (S — 0) case, and only 8.2 m/s for the optimised case, so that the ratio 
is about 5.6, which seems interesting to note. The relative velocity is shown in Figure where its peak to peak value 
is apparent. It is seen from the figure that there are long periods of constant Doppler shifts for the optimised case. 
These results could be useful in the context of TDI. 



-5 
-10 
-15 - 
-20 

-25 



5,=0 
5, =5/8 




0.1 



0.2 



0.3 



0.4 0.5 

Time in years 



0.6 



0.7 



0.8 



0.9 



FIG. 5: The relative velocity ^(^12) is shown for the two cases Si = (dashed curve) and the optimised case Si = 5/8 (solid 
curve) over the course of a year. For the optimised case the gain in peak to peak variation in relative velocity is ~ 5.6. 

III. THE EQUATIONS OF MOTION OF SPACECRAFT TO OCTUPOLE ORDER 

A. Geodesic deviation equation to second order in the separation vector 

We begin with the Newtonian geodesic deviation equation to the first order in the separation vector: 

rj a + K^ b = 0, (24) 

where f] a (t), are the components of the 3- vector separating the two particles at any time t and Kg = The 
<E>(r) = —k/r, k = GM, is the Newtonian gravitational potential, where r is radial coordinate with the centre offeree 
at the origin, G the Newtonian gravitational constant and M the mass of the Sun. Here the commas denote partial 
derivatives with respect to the coordinates. Upper and lower indices have been used to facilitate the summation 
convention. 

If the reference orbit is circular with radius R, the geodesic deviation equation in the barycentric frame is given by: 



(25) 



where X a is the vector (X, Y, Z) . In vector notation (a bar over the quantity denotes a vector quantity and a hat 
denotes a unit vector), 



r/ + Sl 2 f}-3Sl 2 R{R.fj) = 0, 



(26) 



where SI 2 = k/R 3 and R denotes a unit vector in the radial direction. 

As seen before, the CW frame is a suitable frame. In this frame, the separation vector fj is related to the CW 
coordinates by, 



x = r\ x cos Sit + if sin Qt , 
y = —rf sin Sit + r\ v cos Sit , 

Z = 7] Z , 



(27) 



Then transforming Eo. ll2l)l) to CW frame, it is found to be identical to the CW equations: 

x - 2fly - 3il 2 x = , 
y - 2fti = , 

z + ft 2 z = 0. (28) 

However, in order to obtain the flexing of the arms we need to go to higher order in a. We go to the second order 
in a or which is equivalent to the second order in separation vector fj. The second order equation contains octupolar 
terms of the gravitational field which essentially drive the flexing of the arms. The Taylor expansion of the Newtonian 
gravitational potential to the third order leads to the following equation for the separation vector fj: 

r + ^if + ±*%rfif = 0, (29) 



where 

^X c + S a c X b + S bc X a ) + ^ 
Or in the vector form, the above equation can be written as 



$ 'fac = —BK^Xc + S a c X b + S bc X a ) + —X a X b X c . (30) 



■jQ 2 3 i ^ 

fj + Q 2 fj ~ 3r! 2 i?(i?.r7) - —(R.fj)fj ~ + ^ 2 R(R.fj) 2 = 0. (31) 

The octupole terms in the equation are quadratic in the separation vector. These equations take on a simple form 
when transformed to the CW frame described earlier in Eq. 0. 
In the CW frame, Eq. 13111 simplifies to: 
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£-2n y -3{l 2 x + —{2x 2 -y 2 - z 2 ) = 0, 

3ft 2 

y + 2flx —xy = 0, 

3ft 2 

z + ft 2 z — xz = 0. (32) 

R 

These equations can also be thought of as higher (second) order CW equations. 

In order to display the order of the various terms in a, it is useful to employ dimensionless coordinates. We define 
the dimensionless variables in units of I as, 

~x„y„z„p 

x=-,y = -,z = -,p=-, (33) 

where p 2 — x 2 + y 2 + z 2 . 

We also set ft to be unity for convenience. In the dimensionless variables, the equations ifTT^f) reduce to: 

x - 2y - 3i + 3a(3x 2 - p 2 ) = , 
y + 2x — 6axy = , 
z + z — 6axz = . (34) 

B. The Equations of Motion via the Lagrangian 

The equations of motion can also derived with the Lagrangian for a central force expanded around a point on the 
reference orbit. The Lagrangian is: 



\ (X 2 + Y 2 + Y 2 ) + . (35) 



Here the mass of the test particle is ignored; it is not important for the central force problem. However if control 
terms and force terms are added then the mass of the test particle must be included. Here we ignore these effects. 
We carry out our analysis in the CW frame. The Lagrangian given by Eq. flHH m the CW-frame takes the form: 



£ = - \± 2 + y 2 + z 2 + (R + xf + y 2 \ + y (R, + x) - iy + 



k 



(R + x) +y 2 + z 2 



1/2 



(36) 



Then the force free Lagrangian equations of motion are given by: 



x-2y- (1- (3){R + x) = 0, 
y + 2x-(l-(3)y = 0, 
'z + (3z = 0, 



(37) 



where, 



= 



2x fx 

1+ r + (r 



R 



Expanded to the octupole order, the quantity is given by: 



'Sx 
~R 



12 



(I) 2 - 3 ( 



(I)' 



-3 



(38) 



(39) 



To this order of approximation, Eq. I|37J1 are identical to Eq. i!32|) . 



IV. THE SOLUTIONS TO THE EQUATIONS OF MOTION 

We first obtain a general solution to the equations of motion in terms of arbitrary constants and then describe the 
family of solutions with tilt S ~ o(a), which matches with the exact Keplerian orbits Eq. (2J to second order in a. 
Finally, we show that the variation in arm-length is optimal in two ways: (i) r.m.s sense and (ii) peak to peak sense. 
This is achieved by optimally choosing the arbitrary constants in the general solution. These results are valid in the 
neighbourhood of the exact orbits considered in section II A. 



A. The General Perturbative Solution around 5 = 

Accordingly, we write, x = xq + axi, y = ya +otyi, z = Sq + az\, where, Xq, yo and zq are solutions to the equations 
of motion to quadrupole order. These solutions satisfy certain criteria which we describe below. 

To the first order in a the xo,yo: z o satisfy the CW equations l|28|l . As shown in paper I, the solutions of these 
equations satisfying no drift, no offset and which have the test-particle equidistant from the origin are the following: 

Xq = ip cos(t + O ) , 
y = -p sin(t + O ) , 
V3 

Z = —p COs(t + (j) ), (40) 

where po and <f>o are fixed for a specific orbit. For the LISA spacecraft, po — -j=. The orbits of the spacecraft are 
easily obtained by the appropriate choice of 4>q. For convenience in solving the equations, we choose <j>o to be zero. 
The orbits with non-zero 4>q are easily obtained by translation. 



To the second order in a, we get the equations for Xi,yi and z\ as, 

x\ - 2yi - 3ii + 3('3xl - p 2 ,) = , 
y"i + 2xi - 6x Q y = , 

zi + zi — 6x zo = 0. (41) 

We first solve the y-equation. We substitute the solutions for i and y~o from the CW solutions ijlHjl . so that the 
y-equation becomes: 

3 

yi + 2xi + -plsin2t = 0. (42) 

Integrating we have, 

3 

yi + 2x x - -pi cos 2t = A, (43) 

where A is an integration constant. We substitute y~i in the equation for xi to obtain: 

1 5 

Xi + xi = - cos 2 t + 2A + - . (44) 



Then integrating we get, 



1 5 

xi = B cost + C sin t cos2t + 2A+-. (45) 

24 8 



Finally, integrating the equation for Zi ,we get, 



1 \/3 

zi = E cost + F siat ^cos2tH . (46) 

4a/3 4 v ' 

After integrating the yi -equation, namely, Eq. lH3^l , we collect the solutions: 

1 5 
Xi — B cost + Csint —— cos 2t + 2 A + — , 

24 8 

1 5 

yi = -2Bsmt + 2Ccost+ -sin2t- 3At 1 + D , 

6 4 

1 /3 

Zi = Ecost + Fsint ^cos2tH . (47) 

4\/3 4 v 7 

Here, A, £?, C, D, E and F are the integration constants and i|17jl is the general solution for the perturbations to 
second order in a. These constants can be determined from initial conditions. However, the term linear in t in the 
solution for yi represents an unbounded drift. Since such solutions imply instability, the drift must be removed. This 
is achieved if we set A = — y| . 



B. The family of solutions with varying tilt near 5 = 

In this subsection we show that the family of solutions with varying tilt 6 ~ o(a) is consistent with the general 
solution obtained in Eq. I|47J1 . For this purpose we consider the exact orbits of the spacecraft given by Eq. © - 101. 
We are therefore choosing a specific set of orbits for which the spacecraft distances are constant and equal to I to the 
first order in the eccentricity (or equivalently a). If these orbits are expanded to second order in a (note 5 = aSi 
where Si ~ o(l)), the following equations can be deduced for spacecraft 1 in the CW frame in which the general 



solution Eq. H7jl was obtained: 



— cost + 2Ra 2 

-2^ sin< + 2i?a 2 
V3 



1 (Si 



1 



— — — cos t — — cos It — — 



Si 



24 



sin t H — sin 2t 
6 



24 



z = i?a cos i + 2i?a 2 



1 /c- 1 „ V3 
— = (oi — llcosi = cos2£ H — 

2V3 4^3 4 



(48) 



They match with Eq. 1^3 if the constants are chosen as follows: A = — ^, 5 = j — C = 0, D = 0, E — ^7§(£l — 1) 
and F = 0. 

Constructing the full solution (adding first and second order solutions) from Eq. ffiljl and Eq. lUTI) . we obtain 
identical solutions, 



I I 2 



1 Si , 

cos t cos 2t 



1 



24 



Si I sini H — sin 2i 

6 



24 



1 



1 ; 2 

z \ = t: cos t + — — — — 

2 2R 2^3 



(5i — 1) cost 



1 „ \/3 
— — cos 2t H — 

4^3 4 



(49) 



where we have now switched back to coordinates having dimensions of length and afi, y±, z\ being the coordinates of 
spacecraft 1. The coordinates of spacecraft 2 and 3 in the CW frame are obtained by replacing t by t — 2ir/3 and 
t — Air/ 3 respectively. 



C. Establishing the optimality of the solution 



We have shown that the second order CW equations are sufficient for an accurate description of the breathing of 
LISA. We are therefore allowed to use the general solution found in section IV A for the purpose of optimisation. If 
we compute the distance between SC 1 and SC 2 using the general solution, with the same arbitrary constants for all 
SC, we find: 



1X2 = I 



1+ oa/3 fh 5B + ^ E ^_ jL cos + L W -V3E)cos29 



-I(3C-x/3>)sin20-^ 
2 48 



The time varying part of l±% is now: 



1X2 -(hi) = 



-—cos8 + -(3B-V3E) cos 2(9+ -<3C - VSF) sin 20- 
16 2 y J 2 K ' 48 



(50) 



(51) 



Since the terms are orthogonal, the minimisation of the variance amounts to minimising each coefficient separately. 
The arbitrary constants, namely, B, C, E, F occur only in the coefficients of the cos 29 and sin 29 terms while the rest 
of the terms have constant coefficients. Thus the minimisation amounts to setting the coefficients of cos 29 and sin 29 
equal to zero. This gives the minimum r.m.s. variation of about 16,000 km which was the same result obtained at 
the end of section II. 

If we now introduce different constants B,C,E,F for each SC, it is easy to show that only terms in cos k6, k = 1,2,3 
and sin 2$ appear in the distance, in such a way that only the coefficients of cos 29 and sin 2$ can be made to vanish. 



The result is then exactly as the preceding one, showing that this is indeed a true r.m.s. optimum in the neighborhood 
of the exact solution that we have assumed. This is a general result. 

We also minimise the peak to peak flexing by considering the general solution. First, we consider the particular 
case when 3C — y/3F — 0. Then as in section II, 



aly/3 

Al 12 = 



x 3 



12(35 - V3E)x 2 + 3x- 12b] , (52) 



24 

where x = cos 9. Following the analysis analogous to section II, we set the derivative of AZ12 with respect to x equal 
to zero which yields the quadratic equation: 

x 2 - 8(35 - ^/ZE)x + 1 = . (53) 

This equation has complex roots when —0.25 < 3B — ^/3E < 0.25 implying that the extrema of Ali 2 are attained 
when x = ±1 and when this condition is satisfied. From Eq. 15211 one obtains the minimum peak to peak variation 
as al/^/3 ~ 48000 km identical to the one found in section II. However, this family is larger, having two parameters 
B and E than the one considered with varying tilt which has just one parameter 5. In particular, we verify that for 
the family of solutions with varying tilt 8 investigated in section IVB, the condition —0.25 < 3B — ^/3E < 0.25, after 
substituting their values in terms of <5i, leads to, 0.5 < <5i < 0.75 in agreement with the result obtained earlier in 
section II. 

If we do not apriori assume that 3C-V3F is zero, then the simple analysis as performed above, fails. We therefore 
resort to numerical methods. We first observe that the time varying part of l\ 2 upto the constant la^/3/A from Eq. 
is given by: 

5 1 

F(a, b;9)=a cos 26 + b sin 29 - - cos 6 cos 39 , (54) 

8 24 

where, a = 3B — ^/3E and b — 3C — %/3F. The peak to peak value of F we denote by AF. We observe that F 
satisfies the following symmetries: F(—a, b; 9) — —F(a, b;n — 9) and F(a, —b; 9) = F(a, b; —9) implying AF(a, —b) — 
AF(—a, b) — AF(a, b). Thus AF is an even function of a, b. We can therefore restrict ourselves to non-negative values 
of a and b when searching for the optimum value of AF. Also for large a, b, AF ~ 2\J a 2 + b 2 , so that the minimum 
of AF must lie in the neighborhood of the origin with a, b ~ o(l). A straight forward numerical computation shows 
that the minimum value of A_F occurs when —0.25 < a < 0.25 and b = 0, giving the minimum peak variation of 
~ 48,000 km as before. This is again a general result. 



V. CONCLUDING REMARKS 

In this work we have studied the variation of the lengths of the arms of LISA in the neighbourhood of the exact 
solution that gives constant armlengths to first order in a. We investigate the flexing of arms of LISA by going to 
the second order in the parameter a. We numerically establish that the second order results provide a near accurate 
description of the orbits and the flexing of the arms. Because the second order expressions are tractable and amenable 
to analytic treatment, a lot of information can be gleaned from the expressions themselves. We show that by slightly 
changing the tilt angle of the plane of LISA from the proposed one of 60° with the ecliptic, one can reduce the flexing 
of the arms. This is investigated with the aid of the exact Keplerian orbits as well as by analysing to second order 
in a. We then obtain the tilt angle for which the flexing of the arms is minimum. We find that the peak to peak 
amplitude of flexing is reduced by a factor of about 2.4. More importantly, in the context of application of TDI 
techniques, we show that the peak to peak Doppler shift is reduced by a factor of 5.6. 

Secondly, we obtain the equations of motion to second order in a from the Newtonian geodesic deviation equation 
to the quadratic order in the separation vector, which involves the octupolar expansion of the Newonian gravitational 



potential of the Sun. These equations can also be obtained with a Lagrangian approach. The general solutions to 
these equations are found perturbatively, in the neighbourhood of the first order solution for which the armlengths are 
constant to first order in a. We then check that the family of solutions with varying tilt is reproduced from the general 
solution by choosing the arbitrary constants appropriately. Finally we use the general solution to establish that our 
solution of the minimum flexing of the arms is more generally valid in the neighborhood of the exact solution assumed. 
This is important from the point of view of a realistic launch because the flexibility in the choice of parameters allows 
for some error, which is inevitable, and this is possible without sacrificing performance. 
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